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Abstract. One interesting class of gravitational radiation sources includes rapidly 
rotating astrophysical objects that encounter dynamical instabilities. We have carried 
out a set of simulations of rotationally induced instabilities in differentially rotating 
polytropes. An n=1.5 polytrope with the Maclaurin rotation law will encounter the 
m=2 bar instability at T/|W| ^ 0.27. Our results indicate that the remnant of this in- 
stability is a persistent bar-like structure that emits a long-lived gravitational radiation 
signal. Furthermore, dynamical instability is shown to occur in n=3.33 polytropes with 
the j-constant rotation law at T/|W| > 0.14. In this case, the dominant mode of insta- 
bility is to=1. Such instability may allow a centrifugally-hung core to begin collapsing 
to neutron star densities on a dynamical timescale. If it occurs in a supermassive star, 
it may produce gravitational radiation detectable by LISA. 

INTRODUCTION 

One interesting class of gravitational radiation sources includes rapidly rotating 
astrophysical objects that encounter dynamical instabilities. Linear stability anal- 
ysis has shown that rapidly rotating bodies may experience global deformations 
due to the growth of unstable azimuthal modes e ±im ^ [1,2]. The mode numbers m 
describe the shape of the induced deformation. For example, an m=l mode could 
result in the development of a one-armed spiral or a simple translation; an m=2 
mode may produce a bar-shaped distortion; an m=3 mode, a triangular distortion; 
and so on. The onset of such an instability occurs when an object's ratio of rota- 
tional kinetic energy T to gravitational potential energy W, (3 = T/\W\, exceeds a 
critical value /3 cr u- 

Both secular and dynamical varieties of these instabilities may exist. A dynamical 
instability is driven by gravitational and hydrodynamical forces and develops on 
the order of the rotation period of the system. A secular instability is driven by 
a dissipative mechanism, such as viscosity or gravitational radiation, and develops 
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on the timescale of that mechanism (which can be many rotation periods). In 
this paper, we focus on the numerical study of dynamical instabilities, since their 
development can be followed in a reasonable amount of computational time with 
explicit hydrodynamical simulations. 

There are several types of astrophysical objects that may encounter these ro- 
tational instabilities. A star that accretes matter and angular momentum from a 
binary companion may be spun up to rapid rotation [3] . A second example is a cen- 
trifugally hung stellar core or "fizzler" [4-6] . The formation of a fizzler begins when 
the core of a massive star depletes its nuclear fuel and begins to collapse to neutron 
star densities. If the core was rotating initially, conservation of angular momentum 
requires that the core spin up as it collapses. This spin up could actually halt the 
collapse if the centrifugal force increases to the point where it overcomes the inward 
gravitational push. The results would be a rapidly rotating, partially collapsed stel- 
lar core. Another example is a cooling supermassive star (mass > 10 6 Mq) that 
also spins up as it contracts [7,8]. Finally, if the merger of a compact binary does 
not result in an immediate collapse to a black hole, the remnant will be a rapidly 
rotating compact object [9]. The nonaxisymmetric deformations induced by rota- 
tional instabilities could result in relatively strong gravitational radiation emission 
from these rapidly rotating objects. 

We have performed Newtonian hydrodynamics simulations of dynamical rota- 
tional instabilities in two different types of objects [10,11]. The much studied m=2 
bar mode is the strongest of the set of global instabilities encountered by an n—1.5 
polytrope (for which the equation of state gives the pressure P in terms of the den- 
sity p as P = Kp( 1+1 l n \ where K is the polytropic constant) with the Maclaurin 
differential rotation law. This instability sets in when (3 ^ 0.27. Our results indi- 
cate that dynamical instabilities also occur in objects with much lower values of f3. 
In fact, n=3.33 polytropes with the j-constant rotation law encounter a dominant 
m=l instability when (3 reaches 0.14. Simulations of the bar mode instability and 
the m=l instability will be discussed in the following sections. 

THE BAR MODE INSTABILITY 

There has been a discrepancy in the outcome of previous hydrodynamical simu- 
lations of the dynamical bar mode. Namely, simulations carried out by Centrella's 
group showed that an object deformed by the bar mode would become axisym- 
metric again after a short interval [12,13]. This is in contrast to the simulations 
performed by several other groups, which resulted in bar-shaped final configurations 
(e.g., [14-16]). 

The degree of asymmetry in the final configuration is important because an 
axisymmetric object (rotating about its short symmetry axis) will not emit gravi- 
tational radiation. If the simulations performed by Centrella's groups's simulations 
are correct, the instability would produce only a short burst of radiation. However, 
if the object retains a bar-like structure, as the other simulations seem to indicate, 
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FIGURE 1. The h + polarization of the gravitational waveform from a V code bar mode simu- 
lation, as viewed by an observer located along the rotation axis at a distance r from the system. 
The value rh + has been normalized to R~ 1 c~ A M 2 G 2 \ the time is normalized to the dynamical 
time t D = [R 3 /(GM)] 1 / 2 . 
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FIGURE 2. Same as Figure 1, for a simulation performed with the C code. 

it will go on emitting gravitational radiation, producing a longer-lived signal that 
would be easier to detect. 

This discrepancy is illustrated by the gravitational waveforms shown in Figures 
1 and 2. The waveform in Figure 1 is from a simulation performed by Centrella's 
group [12], with a finite-difference hydrodynamics (FDH) code developed at Drexel 
University (hereafter called the V code). This waveform is burstlike; the amplitude 
dies off as the object loses its bar-like shape. The waveform in Figure 2 is from a 
simulation performed by New [14] , with a FDH code developed at Louisiana State 
University (hereafter called the £ code). This waveform rings for the duration of 
the simulation because the object retains its bar-shaped structure. 

Again, the resolution of this discrepancy is important because if the gravitational 
radiation signal persists, the total accumulated amplitude could make it easier to 
detect. The characteristics of New's waveform are such that if the star's initial 
mass and radius are M = IAMq and R = 36 km, and it is located at a distance 
of 20 Mpc, the maximum amplitude is h max ~ 10~ 23 and the frequency f gw ranges 



from ~ 600 — 670 Hz. The scaling is such that h max oc R 1 and f gw oc R 3 / 2 . This 
signal could possibly be detected with advanced interferometers like LIGO II. 

In a recent paper [11], we demonstrate that a good deal of progress has been 
made in resolving the discrepancy among the outcomes of the previous simulations 
of the bar instability. In an attempt to get at the root of the problem, we examined 
the differences between the simulations performed with the V and C codes. One 
of the main differences is that New's simulation with the C code [14] used a con- 
dition called 7r-symmetry. The 7r-symmetry condition enforces periodic azimuthal 
symmetry and thus allows the azimuthal resolution to be doubled as the solution 
only needs to be evolved from to ir. This symmetry condition does not allow the 
growth of odd modes in a simulation. However, the use of 7r-symmetry in a bar 
mode simulation seems justified by analytic perturbation anaylsis, which indicates 
that odd modes will not grow if j3 < 0.32 [17] (/3=0.30 in the initial models used 
by [12-14]). 

In order to investigate if the use of 7r-symmetry was responsible for differences 
in the outcomes of the V and C codes' bar mode simulations, we reran the C 
code simulation with the 7r-symmetry condition turned off (hereafter referred to 
as simulation CI). The CI simulation started with the same axisymmetric initial 
model used in the previous simulations [12,14]. It is constructed in hydrostatic 
equilibrium with Hachisu's Self-Consistent Field (SCF) method [18]. It has an 
n=1.5 polytropic equation of state and is differentially rotating, with an angular 
momentum distribution identical to that of a Maclaurin spheroid [19]. It is a highly 
flattened, rapidly rotating object with (3 ~ 0.3 (recall f3 cr n ~ 0.27). 

Equatorial density contours from the last half of the CI simulation are shown 
in Figure 3. As the final frame indicates, the object appears rather symmetric at 
the end of the run. The resulting waveform (Figure 4) reflects the object's return 
to symmetry. This waveform looks much more like the burstlike waveforms of 
Centrella's group [12,13] (see, e.g., Figure 1). 

Thus the CI run appears to be in better agreement with the V code simulations. 
But the important question is whether the result on which they agree is the "phys- 
ical" result. That is, can we with confidence tell the gravitational wave detection 
community that a star that encounters the bar instability will emit only a burst of 
radiation? 

The answer to that question is no. This is because of an effect that appears late 
in the CI simulation. At about the time that the object loses its bar-like structure, 
the center of mass of the object moves off the center of the grid. This can be 
seen in the final frame of Figure 3. In the absence of external forces, the center 
of mass should not move. This is thus a nonphysical result and must be due to a 
shortcoming in the numerics of the C code. 

Before we could say for certain that an axisymmetric configuration is the correct 
remnant of the bar instability, we needed to find a way to minimize the center 
of mass motion and see how this affected the outcome of the evolution. Further 
testing indicated that increasing the radial resolution of the simulation does indeed 




FIGURE 3. Density contours in the equatorial plane for the later stages of the CI simulation 
are shown. The maximum density has been normalized to unity at t—0 and the contour levels 
are at 0.5, 0.05, 0.005, and 0.0005. 
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FIGURE 4. The gravitational waveform from the CI simulation is shown. The normalizations 
are identical to those used in Figure 1. 



delay the onset of the motion of the center of mass. 

We ran the bar mode simulation with the C code once again, this time with double 
the radial (and axial) resolution used in the CI run (this high resolution run will 
be referred to as C2). In the C2 evolution, the bar-like structure is maintained for 
about 15 more dynamical times (where to — (R 3 /GM) 1 / 2 ) than in the CI run. The 
waveform from this run, shown in Figure 6, has a correspondingly longer duration. 
Doubling the resolution delayed the onset of the center of mass motion, but it did 
not prevent it. Equatorial density contour plots from the latter portion of the C2 
run are shown in Figure 5. As is evident in the final frame, the bar-like structure 
(and waveform) decay at the time that the center of mass starts to move. 

We have demonstrated that the longer the center of mass motion is suppressed, 
the longer the bar-like structure is maintained. Since the center of mass motion 
is unphysical, the decay of the bar is as well. Thus a simulation with very high 
resolution, or better finite-differencing algorithms that did not develop center of 
mass motion, should produce a long-lived nonaxisymmetric structure. A recent 
paper by Brown confirms this [20]. The FDH PPM hydrocode Brown used to 
simulate the bar instability explicitly conserves linear momentum, thus preventing 
center of mass motion from developing. The outcome of Brown's simulation was 
indeed a persistent bar. 

Recall that this is the outcome of the C code's simulation with 7r-symmetry [14]. 
This symmetry condition prevents the growth of any center of mass motion because 
it will not allow an m—1 mode to grow. Our results and the recent paper by Brown 
[11,20] indicate that the physically accurate outcome of the bar instability in this 
object is a persistent bar-like configuration that emits gravitational radiation over 
many cycles and thus may be capable of producing a detectable signal. 
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FIGURE 5. Density contours in the equatorial plane for the later stages of the £2 simulation 
are shown. The contour levels are the same as those shown in Figure 3. 
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FIGURE 6. The gravitational waveform from the £2 simulation is shown. The normalizations 
are identical to those used in Figure 1. 



We have recently also investigated dynamical rotational instabilities in objects 
with lower values of (5 [10]. This study was motivated by the fact that centrifugally 
hung stellar cores, or fizzlers, are likely to have f3 < 0.27 and thus will not encounter 
the particular bar mode instability that was discussed in the preceeding section 



Previous authors have identified dynamical instabilities in toroidal configurations 
at values of (3 ^ 0.1 [23,24]. Thus we decided to investigate the stability properties 
of spheroidal configurations with off-center density maxima (thinking they may be- 
have more like tori). Our goal was to determine whether spheroidal configurations, 
representative of fizzlers, could encounter instabilities at lower values of (5. 

The initial models for this stability study were constructed with Hachisu's SCF 
method. An n=3.33 polytropic equation of state was used and is representative of a 
partially collapsed stellar core. Because we wanted to study models with off-center 
density maxima, we needed to use a rotation law other than the Maclaurin law. 
The j-constant law does produce off-center density maxima in n=3.33 polytropes. 
The angular velocity distribution Q(zu) given by this law is 



where w is the cylindrical radius and d is an arbitrary constant. As d approaches 
zero, the specific angular momentum approaches the constant jo- We constructed 
models with d=0.2. Figure 7 displays meridional density contour plots of four of 
the equilibrium models we constructed. Off-center density maxima appear as (5 is 
increased. 

We performed simulations of these models with the C code and with Brown's 
PPM hydrocode. The £ code uses a cylindrical grid and Brown's code uses a 
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FIGURE 7. Density contours in the meridional plane are shown for 4 equilibrium models with 
d=0.2 and n=3.33. The models are labeled with their values of (i = T/\W\. The initial maximum 
density has been normalized to unity and the contour levels are at 0.9, 0.1, 0.01, and 0.001. 



Cartesian grid. We have evolved all four of the models shown in Figure 7. The 
onset of instability is near /3=0. 14. Both the /3=0.14 and 0.18 models are clearly 
unstable (see below). We have also run simulations of the models with (3=0.12 and 
0.09. The (3=0.12 model was run for about 40 t D . At the end of the simulations, 
the m=l mode was just starting to grow. The /3=0.09 model was run for 35 tu and 
showed no mode growth. We plan to run longer, higher resolution simulations to 
determine the value of (3 at which the instability sets in. Note that no significant 
center of mass motion was seen in these simulations. 

The evolution of the model with (3=0. 1A exhibits a dynamical instability that is 
very similar in simulations performed with both codes. The development of the 
instability is shown in the equatorial density contour plots of Figure 8 (from the C 
code run). The torus pinches off to produce a single high density region, indicative 
of a m=l mode. This dense region starts to collapse at late times, as is evident from 
the appearance of higher density contours. Figures 9 and 10 show the amplitude of 
the modes that grow during the evolutions performed with both codes. The m=2, 
3, and 4 modes grow in sequence after the m=l mode. The pattern speeds for the 
m=l and m=2 modes are identical. Thus the m=2 mode is an harmonic of the 
m=l mode. The constant amplitude of the m=4 mode in Brown's simulation is a 
result of the symmetry of his Cartesian grid. 

The evolutions of the (5=0.18 model also exhibit dynamical instability. The mode 
amplitude plots from these runs are shown in Figures 11 and 12. The instability 



FIGURE 8. Density contours in the equatorial plane are shown for the model with /3=0.14. 
The darkness of the shading increases as the density decreases. The contour levels are 0.01, 0.1, 
0.9, 2, and 4 times the maximum density at the initial time. 
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FIGURE 9. The growth of the amplitudes \A m \ for m=l (thick solid line), m=2 (thin solid 
line), m=3 (dashed line), and m=4 (dotted line) is shown for the C code simulation of the /3=0.14 
model. These amplitudes were calculated in the equatorial plane for a ring with radius m=0.32. 




FIGURE 10. Same as Figure 9, for Brown's simulation. 
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FIGURE 11. Same as Figure 9, except that the model has /3=0.18. 
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FIGURE 12. Same as Figure 10, except that the model has /3=0.18. 

grows more quickly in this model than in the (3=0.14 model, as expected. The m=l 
mode grows at about the same rate in the simulations of both codes. The C code 
run shows the m=2, 3, and 4 modes growing in sequence, just like in the /9=0.14 
simulation; the m=2 mode is again an harmonic of the m=l mode. Brown's run, 
however, shows the m=l and m=2 modes growing at the same rate. These modes 
are not harmonics in his simulation. We plan to investigate the differences in these 
results with higher resolution simulations. 

Our study demonstrates that dynamical instabilities can occur in differentially 
rotating polytropes with relatively low values of f3. Note that Pickett, Durisen, 
and Davis [15] also found a dominant m=l instability in a centrally condensed 
configuration. The instability in their model set in at /3=0.2. The model had 
an ra=1.5 polytropic equation of state and was constructed with the so-called n'=2 
rotation law. This rotation law places more angular momentum in the outer regions 
of the model than does the j-constant law. Their model did not have off-center 
density maxima. 

We plan to carry out more detailed studies to investigate the character of these 




unstable modes and their properties for various values of the polytropic index n 
and the rotation law parameter d. 

If this instability occurs in centrifugally hung stellar cores, collapse to neutron 
star densities may result. Our simulations do show that the density is increasing at 
the end of the unstable runs. Further studies are needed to determine how dense 
the remnant becomes and to see if the m—1 mode will result in the star moving 
with velocities comparable to those of actual neutron stars. 

Longer runs on larger grids are needed to obtain the full gravitational radiation 
waveforms. However, we can estimate the properties of the emission during the 
initial stages of the development of the instability. For a fizzler that starts out with 
M ~ 1.4M© and R ~ 200km, the peak emission will occur at f gw ~ 200Hz. The 
peak amplitude h max will be ~ lO -24 ^ 1 for (3 = 0.14, and ~ lO -23 ^ 1 for (3 = 0.18. 
Here, r 2 o is the distance to the source in units of 20 Mpc. Emission from such 
unstable cores may be detectable with advanced ground-based interferometers like 
LIGO II. 

An even more optimistic scenario for the detection of gravitational radiation 
occurs if this type of instability is encountered by a cooling and contracting super- 
massive star [7,8]. If the star's ratio GM/Rc 2 is ~ 15 when the instability develops 
(a value that is approximately appropriate for an uniformly rotating star), f gw 
will be ~ 10~ 3 Hz and h max will be ~ lO" 18 ^ 1 for ^ = 0.14 and ~ lO" 17 ^ 1 
for (3 = 0.18. Such signals would be easily detectable by the space-based LISA 
detector. 
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